Simulated three-component granular segregation in a rotating drum 
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Discrete particle simulations are used to model segregation in granular mixtures of three different 
particle species in a horizontal rotating drum. Axial band formation is observed, with medium-size 
particles tending to be located between alternating bands of big and small particles. Partial radial 
segregation also appears; it precedes the axial segregation and is characterized by an inner core 
region richer in small particles. Axial bands are seen to merge during the long simulation runs, 
leading to a coarsening of the band pattern; the relocation of particles involved in one such merging 
event is examined. Overall, the behavior is similar to experiment and represents a generalization of 
what occurs in the simpler two-component mixture. 

PACS numbers: 45.70.Mg, 45.70.Qj, 64.75.+g, 02.70.Ns 



I. INTRODUCTION 

One of the more fascinating properties of granular mat- 
ter is the ability of appropriately driven mixtures to sepa- 
rate into their individual components, despite the appar- 
ent lack of energetic or entropic advantages of an unmixed 
state. The segregation of a binary mixture contained in 
a partially filled, horizontal, rotating drum is an exten- 
sively studied problem of this class, one with obvious 
industrial importance. The most prominent characteris- 
tic exhibited by this system is axial segregation, in which 
a pattern of bands of alternating particle species forms 
along the cylinder axis; the bands merge and the pat- 
tern coarsens over time, at a rate that drops to such a 
low level that it is unknown whether the state eventually 
reached is stable, or merely longlived. A second form of 
segregation, though one that can require greater effort to 
observe, occurs in the radial direction, producing a core 
rich in small particles surrounded by an outer layer of 
mainly big particles; this core can either be a transient 
on the path to later axial segregation, or a feature that 
persists even when full axial segregation is apparent from 
external views. 

In the case of mixtures of two species, early experimen- 
tal efforts involved direct observation of the axial band 
structure at the outer surface [l|, whereas some of the 
subsequent studies succeeded in examining the interior 
using MRI 0, [1[ and relating the development of the 
axial bands to bulges occurring in the radial core. Time- 
dependent behavior can occur subsequent to the initial 
appearance of the axial bands; narrower bands gradually 
merge to form broader bands [J| and traveling surface 
waves have been noted [f|. Complicating the behavior 
are results showing that the ratio of cylinder to parti- 
cle diameter determines if axial segregation can occur 
and whether its appearance depends reversibly on rota- 
tion rate [f| . Further recent experimental examples 0, HI 
reveal more about the richness of the segregation effect 
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and the associated dynamics under wet conditions; these 
are relevant because the behavior of slurries tends to be 
similar to dry mixtures. Overall, no consensus has yet 
emerged as to the nature of the mechanisms underlying 
this complex set of segregation phenomena. For further 
discussion of the difficulties inherent in understanding 
granular matter in general see, e.g., Refs. [9l. [ToL Hl| . 

In contrast to the many studies involving two- 
component systems (only a selection of which were cited 
above), three-component systems have received far less 
attention. The experimental results available for systems 
with three particle species [l2| were obtained using glass 
beads of different sizes (an earlier brief discussion appears 
in [Hj]). The description of the experiment relates that, 
after beginning with an initially uniform mixture, the 
first effect to occur is the appearance of radial segrega- 
tion, in which there is an inner core of small particles and 
an outer layer of big particles, with the medium-size par- 
ticles lying predominantly in between. This is followed by 
the development of axial bands of medium-size particles 
in the mainly big-particle outer region, followed by the 
appearance of small-particle bands within the medium- 
particle bands. Replacing the opaque big and/or medium 
particles with transparent equivalents reveals that the 
small-particle core persists after the onset of axial seg- 
regation, and that a layer of medium-size particles sur- 
rounds it. Systems with more than three components are 
also mentioned, with four components showing behavior 
analogous to three, while beyond four the axial segre- 
gation is no longer seen, but these lie outside the scope 
of the present work. As with two-component systems, 
the behavior is certain to be influenced by the relative 
species sizes and concentrations, as well as by the many 
parameters required to fully specify the system, includ- 
ing the fill level, frictional properties and rotation speed, 
although their influence on the nature of the segregation 
is not described. 

Simulations of granular systems based on discrete- 
particle models have proved capable of reproducing both 
the axial and radial forms of segregation, in qualitative 
agreement with experiment. The present paper extends 
an earlier study [14| . in which two-component systems 
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were modeled. In that paper it was shown how the ap- 
pearance of radial segregation preceded axial segregation, 
and that depending on the circumstances, a radial core 
of small particles could still be present even after the 
exterior view showed essentially complete axial segrega- 
tion. Owing to the need for long runs and large systems, 
only a limited set of parameter combinations could be 
considered. The problem of following the development 
of segregation in detail is further exacerbated by a lack 
of reproducibility between runs differing only in the ran- 
dom velocities assigned to the particles in the initial state 
that result in different segregation scenarios. While sim- 
ilar effects occur experimentally, the implication is that 
multiple simulation runs are essential to ensure that the 
'typical' behavior is captured. In view of the potentially 
heavy computational demands, a thorough study of the 
two-component system and a detailed comparison with 
the variety of experimental data available has yet to be 
undertaken. 

The obvious next stage in this exploration is to increase 
the number of granular species, and while this has been 
carried out experimentally, as indicated above, the cor- 
responding simulations have yet to be carried out. The 
goal of the present paper, therefore, is to demonstrate 
that the model used previously is also able to produce the 
appropriate kinds of segregation even when three species 
are involved. Given the very limited amount of experi- 
mental information presently available for this problem, 
any systematic coverage of the parameter space is prema- 
ture; the discussion is therefore confined to specific ex- 
amples illustrating the main features of simulated three- 
component segregation. 

II. METHODOLOGY 

The granular model employs spherical particles whose 
overlap repulsion is based on a continuous potential [H, 
EH E3 1 exactly as in [3] • The normal force between a 
pair of particles includes a repulsion that depends 
linearly on the overlap and a velocity-dependent damping 
force, 

fn [kn{d%j Vij i) ^ni^ij ' ^iji Tij ^ dij ) (-0 

where and tiy are the particle separation and relative 
velocity, and dij = (di + dj)/2 is the mean diameter. 
The two parameters appearing in Eq. |T]), k n governing 
the stiffness and 7„ the normal damping, are the same 
for all particles. 

Sliding of particles that are within interaction range 
is opposed by both frictional damping and a transverse 
restoring force, 

ft ' < K [ «£(r)dr, (2) 

J (coll) 

where u*- is the relative transverse velocity allowing for 
particle rotation, and the integral is evaluated as the sum 



of vector displacements over the collision period (addi- 
tionally, whenever the integrated displacement exceeds 
0.1 it is reset to zero). The magnitudes of both terms in 
Eq. arc bounded by l^ CiCj \f n \- The transverse damp- 
ing and static friction coefficients, 7s° C3 and /i CiCj , depend 
on the particle types Cj and Cj ; the transverse stiffness k g 
is the same for all particles. 

The majority of parameter settings are taken from [Tij 
and, as before, the curved cylinder wall and the flat end 
plates are treated as rough and smooth boundaries, re- 
spectively. The value of the gravitational acceleration, 
g = 5, relates the dimensionless MD (molecular dynam- 
ics) units of the simulation to the corresponding physical 
units. If Lmd is the length unit (in mm), then the time 
unit is Tmd ~ 10 -\/ELmd s. A cylinder rotating with 
angular velocity fl (MD units) has an actual rotation 
rate of ^1/(2ttTmd) ~ 7 .10,/ \/Lmd Hz; for 4mm parti- 
cles, O = 0.2 is equivalent to an experimental 43rpm. 
Other details of the simulations are covered in |18l|. and 
the required MD techniques are described in fl9j ]. 

Efficiently evaluating the displacement sums in Eq. ([2]) 
requires maintaining a list of all currently interacting par- 
ticle pairs in a form that is readily accessible given the 
identity of the pair involved, but also easily modified as 
pairs are added or deleted due to the continually chang- 
ing identities of interacting neighbors. Of the approaches 
available for organizing the data, the one used here (as 
well as in [l4|, although the description was omitted) is 
the following. 

A table is constructed that holds a set of linked lists 
of interacting pairs (i, j), ordered so that i < j. Each en- 
try consists of the value of j, a pointer to another entry 
used to link entries with a common i, the current value of 
Uij — J2(coU) v fj( T ) needed in Eq. ([2]), and the timestep 
number t^ at which iiy was last updated. There is a sep- 
arate list for each particle i, and a pointer hi to the first 
entry in its list; from hi the list can be scanned using the 
pointers to obtain all interacting neighbors with j > i, 
until a null pointer ends the list. For each interacting 
pair encountered during the force calculations, the corre- 
sponding and tij are updated; new entries are readily 
added as necessary, whereas if a previously interacting 
pair is no longer in range its entry will not be updated. 
After all the forces have been evaluated, those entries 
whose are not current (because > dij) have their 
Uij reset to zero. Expired entries need not be deleted im- 
mediately since pairs are capable of moving in and out of 
range many times, in which case, if the relevant entry is 
still in the table it will be reused. (On a shared-memory 
computer, if the neighbor list references the relevant ta- 
ble entries, the force computations can be carried out in 
parallel [19] , even with this history-dependent interac- 
tion.) 

The force parameters are k n = 1000, 7„ = 5, and 
kg = 500. Transverse damping (or dynamic friction) de- 
pends on particle type; if 6, m and s denote big, medium 
and small sizes, then 7^ = 10, = 6, 7f s = 2, and 

the lower of the values is used for unlike species (the bb 
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and ss values are from [14|, and the mm value is their 
mean). Finally, the relative values of the static friction 
coefficients, e.g., p bb /p ss , are set equal to the ratio of the 
corresponding 7^ values, with the larger of the pair equal 
to 0.5. In general, there seems to be considerable latitude 
in choosing the parameter settings, and what is required 
at this exploratory stage is the ability to produce results 
that are qualitatively reasonable. 

The nominal particle diameters correspond to the in- 
teraction cutoff. As in for small particles d s — 
2 1 / 6 ps 1.122 (MD units)'. For big particles d b = 2d s , 
and for medium particles d m = 1.3d a ; this choice of rela- 
tive sizes is not too different from that used experimen- 
tally. The actual small-particle diameters are uniformly 
distributed over a narrow range [d s — 0.2, d s ] (the mean 
diameter of the small particles is then close to unity), 
and likewise for the other species. Particles all have the 
same material density. 



III. RESULTS 
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FIG. 1: (Color online) Axial space-time plot for run #A 
(L = 240, L/D = 8, riR — 4260); red, green and blue (or 
medium, light and dark shades of gray) denote whether the 
big, medium or small particles have the highest volume frac- 
tion. Time is along the horizontal axis, and the vertical axis 
corresponds to axial position. 



Three simulation runs are singled out for analysis here, 
each aimed at demonstrating particular features of the 
behavior. They employ cylinders with diameter D = 30 
or 40 (reduced units) that are rotated with angular ve- 
locity u> = 0.2 where particle flow is continuous; these 
are the values used for the majority of runs in the two- 
component study [HI . Runs are begun with the particles 
arranged on a lattice at a specified density (this deter- 
mines the eventual fill level) and assigned random veloc- 
ities. Particle species is assigned randomly, producing 
an initially mixed state; the relative populations of the 
three kinds of particles are chosen to ensure equal volume 
fractions of all species. 

The first of the runs, #A, involves a cylinder of length 
L = 240 and diameter D — 30 (aspect ratio L/D = 8); 
the initial filling density is p — 0.25, resulting in a sys- 
tem with N — 29 704 particles. The total run length is 
riR : = 4260 revolutions, and the purpose of this compar- 
atively long run (by simulation standards) is to demon- 
strate the ability to produce an apparently stable ax- 
ial band pattern after the effects of the initial transients 
have disappeared. The second run, #5, employs a longer 
cylinder, L = 360, with the same D and filling density, 
so that N — 44 556; the run length ur — 2300 is shorter, 
and the higher aspect ratio (= 12) can accommodate a 
more extensive set of axial bands. In the third run, #=C, 
the cylinder length is reduced to L = 120, the diameter 
increased to D = 40, and the fill level raised by setting 
the initial density to p — 0.4 so that N = 47 520; this is 
a shorter run whose goal is to provide increased space for 
observing the development of radial segregation. 

The behavior encountered in each of these runs is de- 
scribed in detail below. This is accomplished using a 
combination of axial space-time plots that summarize the 
overall development of axial segregation, graphs showing 
the concentrations of each of the particle species along 



the axial and radial directions, and images of individual 
configurations that reveal where particles actually reside. 

Figure Q] shows the axial space-time plot for run j/A. 
In view of the need to distinguish between three species, 
rather than just two, the color coding shows only the 
volume- weighted majority species, irrespective of the de- 
gree to which the other species might be present; in this 
respect the present space-time plots differ from [14| where 
a continuous color gradient designates the relative vol- 
ume fraction itself. The plot shows that during slightly 
over half of the run a series of band merging events oc- 
curs, most within the first 1200 revolutions and the last 
at approximately 2500, eventually producing a state with 
10 axial bands (counting all species, where one of the end 
bands is only weakly defined), and there is no apparent 
change thereafter. While the possibility of future change 
cannot be excluded, none was observed when the run was 
subsequently extended by an extra 1000 revolutions (not 
shown). The rapid, small fluctuations at the band bound- 
aries are artifacts of the discrete nature of the majority 
rule by which color is determined. 

More detailed information concerning the band occu- 
pancy at the end of run #A appears in Fig. [2] The plot 
shows the volume fractions of each of the particle species 
in a sequence of slices oriented normal to the cylinder 
axis; fluctuations in the measurements are reduced by av- 
eraging over ten successive configurations, sampled once 
per revolution. It is clear that the b and s particles are 
well segregated axially. On the other hand, the axial sep- 
aration of 771 and b particles is less pronounced than for m 
and s; while the m particles exhibit a distinct preference 
for the regions located between b and s bands, and their 
occupancy drops to a very low level in the middle of the 
s bands, they are seen to maintain a significant minority 
presence throughout the b bands. 
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FIG. 2: Volume fractions of the three particle species mea- 
sured in slices normal to the axis for the axially segregated 
state near the end of run #A (distance is expressed in dimen- 
sionless MD units). 



FIG. 4: Volume fractions of the particle species early in run 
#A measured in slices parallel to the free surface (the cylinder 
axis is at the origin). 




FIG. 3: Development of the peaks in the the axial distribution 
of small particles for run #A. 

Adding an extra time dimension to Fig. [2] provides a 
more detailed description of the segregation process than 
is available in the space-time plot, but the results for 
the three species must be graphed separately to ensure 
the features are visible. The time-dependent distribution 
of small particles is displayed as a surface plot in Fig. O 
prominent features include the emergence of the peaks in 
the axial distribution and the subsequent peak changes 
corresponding to axial band reorganization. 

Figure 0] shows the volume fractions at an early stage 
of run #A (after approximately 37 revolutions, prior to 
the appearance of any axial effects) in a manner aimed at 
revealing the presence of radial segregation. In order to 
achieve this, the system is divided into a series of slices 
that are oriented parallel to the mean direction of the 
upper free surface of the mixture (as in [lij , the slope is 



determined by a linear fit to the inner 2/3 of the surface 
away from the curved boundary) and the relative volume 
fractions in the slices evaluated (also limited to a region 
of the same width to avoid bias); the results are again 
averaged over ten configurations, at one-revolution inter- 
vals, for improved statistics. The b and s particles are 
seen to experience a certain degree of radial segregation, 
with s dominating the interior and b on the outside, but 
the m particles exhibit no radial preference (the behavior 
is slightly different in run #C, below, where the particle 
layer has double the maximum thickness due to larger D 
and a higher fill level). 

Figure [5] shows the axial space-time plot for run #5, 
where the longer cylinder allows the formation of a 
greater number of axial bands. Several band merging 
events are apparent, most during the first 400 revolu- 
tions and the last at about 1200, resulting in 19 bands at 
the end of the run; the final band merge will be examined 
in greater detail below. A plot of the volume fractions 
as a function of axial position at the end of the run is 
shown in Fig. [SJ once again the b and s particles are well 
segregated, while the m particles exhibit the same pref- 
erence for the regions between the b and s bands and a 
minority presence inside the b bands, as in run #A. 

While the colored space-time plots only allocate a sin- 
gle column of pixels to describe the state of the entire 
system at each instant, images of the evolving patterns 
viewed from above provide a more detailed history resem- 
bling what would be seen experimentally. Figure [7] shows 
the state of run #5 after approximately 20, 60, 100, 200, 
400, 800 and 1600 revolutions; the images cover the early 
portion of the run where there is no hint of axial seg- 
regation (although there are radial processes occurring 
beneath the surface), followed by the pattern coarsening 
stage, and the steady state that persists throughout the 
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FIG. 6: Volume fractions of the particle species in the axially 
segregated state near the end of run =fj=B. 



latter part of the run. 

Figure [8] contains views of the final state of run #5 as 
seen from the outside, along with images of the individ- 
ual species that correspond to what might be observed 
experimentally, either using MRI, or by replacing any 
two of the species with transparent particles having sim- 
ilar mechanical properties. The first two views are of the 
system seen from above and below (the view direction is 
approximately normal to the free surface) ; the remaining 
views each show a single particle species. The b and s 
particles are seen to occupy well-defined axial bands; the 
m particles, however, show only partial axial segregation, 
an observation confirmed quantitatively in Fig. [6j A fea- 
ture apparent from the pictures is that the boundaries 
of the b and s bands are not vertical, with the s bands 
being narrower at the upper free surface. 

The distribution of volume fractions for run af- 
ter 35 revolutions, evaluated as before to demonstrate 
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FIG. 8: (Color online) Final state of run #B; the full system 
seen from above and below, and separate views of the big, 
small and medium size particles, colored copper, gold and 
silver, respectively). 



the presence of radial segregation (also averaged over ten 
configurations), is shown in Fig. [5] Unlike the case of run 
#A earlier, the behavior is more clearly resolved due to 
the larger cylinder diameter and higher fill level. There is 
now also a suggestion of the m particles having preferred 
radial positions, although the effect is still much weaker 
than for the b and s particles. A view showing three slices 
cut out of the system appears in Fig.QIjl what can be seen 
here reinforces the conclusion from Fig. [5J that despite 
the noise, the b and s particles show distinct preferences 
for the exterior and interior, respectively, with the m 
particles also preferentially located towards the exterior, 
though to a lesser degree than b (axial segregation begins 
later in this run, although the low aspect ratio impedes 
its development). 

When comparing the simulation results with experi- 
ment, apart from any possible shortcomings of the model, 
allowance must be made for the relatively small system 
sizes that can affect the results. Given that the compu- 
tational requirements depend on the product of N and 
riR, the need to avoid excessively long simulations means 
that the ratio of cylinder diameter to particle size, D/d s , 
must be smaller than in experiment; in addition to finite- 
size effects in general, this particular limitation has been 
shown to influence behavior experimentally Q. Nev- 
ertheless, to the extent that it is possible to compare 



6 



0.025 



0.020 



0.015 



0.010 



0.005 



0.000 





1 1 1 1 1 1 1 
□ small 




« medium i 






I 


\5 



5 10 
position 



15 



FIG. 9: Volume fractions of each particle species early in run 
# C measured in slices parallel to the free surface (the cylinder 
axis is at the origin). 




FIG. 11: (Color online) Band merging in run H^B\ the views 
show the full system and the individual particles from the b, 
s and m bands involved (only the inner portions of the m 
bands), both before and after the event. 




FIG. 10: (Color online) Run #C after 35 revolutions; three 
narrow slices through the system are shown. 



the three-component simulations with the limited exper- 
imental data available [l2| , the similarity is encouraging, 
in particular, the organization of the axial bands, with 
the bands of m particles located between alternating b 
and s bands. The fact that experiment indicates that a 
core of m particles persists within the b bands is also in 
agreement with the simulations, even if in the latter the 
core is less well defined. 

As always, however, it is the deviations from experi- 
ment that must be considered. For example, experiment 
suggests that, in addition to the axial m bands being 
more sharply defined, there is apparently an innermost 
core of s particles (surrounded by m particles) . Whether 
such differences, in particular the latter, are due to basic 
flaws in the model, or whether they are merely conse- 
quences of a less than ideal parameter choice (perhaps 
as simple as the range of particle sizes) , or too narrow a 
cylinder, must await future study, both simulational and 
experimental. 



The final example considered involves the band merg- 
ing event of run #5 that occurs between revolutions 700 
and 1450 at a location in the cylinder sufficiently far from 
the end plates to be free from any interference. An anal- 
ogous study was carried out for the two-species system 
[14| , but with three species the details are more complex. 
Figure QT] shows what becomes of the occupants of the 
merging b bands and the vanishing s and m bands. The 
first two images are of the entire system seen from above, 
before and after the event; pattern changes are limited 
to the neighborhood of the merging bands. Each subse- 
quent pair of images shows the particles of a single species 
within the affected bands, both before the event (includ- 
ing particles beneath the surface intruding into adjacent 
bands) and afterwards; in the case of the m particles, 
only those initially in the inner halves of the two bands 
are considered. The b bands are seen to merge with only 
the slightest scatter outside the newly formed band. The 
s band disperses into the adjacent bands with practically 
no scatter beyond the bands receiving these particles. In 
the case of the m bands (where only the inner half-bands 
are tracked) the contents are seen to disperse to a greater 
degree than the b bands; such behavior is consistent with 
the fact that axial segregation, and hence confinement, 
is weakest for the m particles. 



IV. CONCLUSION 

Particle-based simulations of three-component granu- 
lar systems in a rotating drum have been used in demon- 
strating that both axial and radial segregation can be 
produced using a simple model for the granular particles 
and their interactions. In the case of the two-component 
systems considered previously, there was already a rel- 
atively large set of parameters whose settings are po- 
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tentially capable of influencing the behavior; with three havior, and developing a closer connection between sim- 
components the set is even larger. Exploring the con- ulation and experiment, are subjects for future study, 
sequences of systematic parameter variations on the be- 
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